function detectPower = evalSeq(transSeq,SOA, HRF_file, HRF_res)
%evalSeq: Calculate the detection power  of a given sequence
%
%
%   Inputs
%   -------
%   transSeq: Sequence of transition distances
%   SOA: Stimulus onset asynchrony, or inter-trial interval (in ms)
%   HRF_file: Path to the hrf .txt file
%   HRF_res: Temporal resolution of the hrf file (in Hz)
%   
%   PS: NA distances in the transition sequence should be represented with -1
%
%   _______________________________________________
%   by Marcelo G Mattar (13/07/2012)
%   mattar@sas.upenn.edu

if nargin<4
    HRF_res = 10;
end

if nargin<3
    HRF_file = 'hrf.txt';
end


rawHRF = dlmread(HRF_file);
HRF_res = (1/HRF_res)*1000; %convert to ms


%% Mean center
% Mean center non-NA values
transSeq_meancentered = transSeq - ((transSeq~=-1) * sum(transSeq(transSeq~=-1))/sum(transSeq~=-1));

% Set NA values to zero
transSeq_meancentered(transSeq==-1) = 0;


%% Upsample transSeq
factor = SOA/HRF_res;
transSeq_upsampled = zeros(length(transSeq_meancentered)*factor,1);
for i=1:length(transSeq_upsampled)
    if mod(i-1,factor)==0
        transSeq_upsampled(i) = transSeq_meancentered((i-1)/factor + 1);
    else
        if ceil((i-1)/factor + 1) < length(transSeq_meancentered)
            transSeq_upsampled(i) = transSeq_meancentered(floor((i-1)/factor + 1)) + ((transSeq_meancentered(ceil((i-1)/factor + 1)) - transSeq_meancentered(floor((i-1)/factor + 1)))/factor) * (mod(i-1,factor));
        end
    end
end

%% Convolve transSeq with the HRF

% FFT both hrf and transSeq
hrf_fft = fft(rawHRF,length(transSeq_upsampled));
transSeq_fft = fft(transSeq_upsampled);

% Normalize hrf_fft
hrf_fft = hrf_fft/max(abs(hrf_fft));

% Apply notch filter
hrf_fft(f_axis(length(hrf_fft),0.1)<0.01) = 0;

% Multiply (convolve) transSeq_fft and hrf_fft
convolved_fft = hrf_fft .* transSeq_fft;

% IFFT the product
convolved = ifft(convolved_fft);

%% Calculate Detection Power
detectPower = var(convolved)/var(transSeq_upsampled);

%==========================================================================





function freq_axis = f_axis(numPoints, sampleRate)
    N = numPoints; 

    % Sampling frequency
    Fs = 1/sampleRate;

    % Calculate the number of unique points
    NumUniquePts = ceil((N+1)/2);

    % Create the first half of the frequency axis
    f_1sthalf = (0:NumUniquePts-1)*Fs/N;

    % Create the second half of the frequency axis (reverse order)
    if rem(N,2)
        f_2ndhalf = f_1sthalf(2:end);
    else
        f_2ndhalf = f_1sthalf(2:(end-1));
    end

    % Attach both halfs together to form the full frequency axis
    freq_axis = [f_1sthalf fliplr(f_2ndhalf)];